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Nuclear quantum effects in water exchange around lithium and fluoride ions 
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We employ classical and ring polymer molecular dynamics simulations to study the effect of nuclear 
quantum fluctuations on the structure and the water exchange dynamics of aqueous solutions of 
lithium and fluoride ions. While we obtain reasonably good agreement with experimental data 
for solutions of lithium by augmenting the Coulombic interactions between the ion and the water 
molecules with a standard Lennard-Jones ion-oxygen potential, the same is not true for solutions of 
fluoride, for which we find that a potential with a softer repulsive wall gives much better agreement. 

A small degree of destabilization of the first hydration shell is found in quantum simulations of both 
ions when compared with classical simulations, with the shell becoming less sharply defined and 
the mean residence time of the water molecules in the shell decreasing. In line with these modest 
differences, we find that the mechanisms of the exchange processes are unaffected by quantization, 
so a classical description of these reactions gives qualitatively correct and quantitatively reasonable 
results. We also find that the quantum effects in solutions of lithium are larger than in solutions 
of fluoride. This is partly due to the stronger interaction of lithium with water molecules, partly 
due to the lighter mass of lithium, and partly due to competing quantum effects in the hydration of 
fluoride, which are absent in the hydration of lithium. 


I. INTRODUCTION 

The exchange of water around solvated ions is a key 
step in bioc hemical processes such as reactions at protein 
surfaces,El^and in geochemical processes such as the depo¬ 
sition of clays and minerals from aqueous solutionJ^ and 
the absorption of contaminants by minerals in soil.^ In 
all three of these examples, understanding the seemingly 
simple reaction 

I*" (H20)„+H20* ^ (H20)„_i (H20*)+H20, (1) 

is a prerequisite for understanding much more complex 
processes. A great deal of work has therefore been done 
on this water exchange reaction. In particular, the mean 
residence time (MRT), r, of water in the hydration shells 
of ions - or equivalently its reciprocal, the exchange rate 
k - has been the subject of numerous experimental and 
computational studies.l^^®^ When the exchange is slow 
enough, the MRT can be inferred from nuclear magnetic 
resonance (NMR) measurements.li^K^When the exchange 
is faster, NMR only provides an upper bound on the 
residence time, which can be reduced to a 100 ps upper 
bound with an incoherent quasi-elastic neutron scattering 
experimentShorter (sub 10 ps) residence times have 
also recently been measured more precisely for certa in 
anions by two-dimensional infrared spectroscopy.^^HHl 

One popular computational route to the MRT was sug¬ 
gested by Impey et alW. t can be calculated as the decay 
time of a residence correlation function. This method, 
often denoted “direct simulation”, relies on being able 
to observe exchange events in a simulation, so cannot 
be used if r is too long. However, for most monova- 
len t and m onatomic ions, r is on the order of 10 ps or 
les^mHH] exception of Li+, for which values 


as high as 400 ps have been reportecPSl) , and the cal¬ 
culation of the MRT in this way is generally viable. A 
similar hydrogen-bo nd cor relation function, introduced by 
Luzar and Chandler ,l2Sl2Sl has also been used to calculate 
reside nce ti mes of water in the first hydration shells of 
anions An alternative method was introduced by Rey 
and Hynes,l2l2ll .^^ho treated the escape of water from the 
first hydration shell of the ion as an activated chemical 
process with rate constant k = calculated using the 
method of reactive flux. This has the advantage that 
there are no restrictions on the magnitude of t, allowing 
the study of exchanges that are too slow for the direct 
simulation method.1121 

There are several other methods that have been used 
to calculate r: a few ab initio studies of water exchange 
have found the MRT by cou nting the number of exchange 
events per unit time .1^21311 Unlike the direct simulation 
method, this does not require the fitting of a correlation 
function to an exponential decay, so may be advantageous 
when this fitting is problematic. Laage and HyneP^ also 
proposed calculating r using the stable-states picture of 
chemical reactions, which is similar to the reactive flux 
method, but can also be used if the MRT does not reach a 
plateau value. Finally, Kerisit and RosscJi^used transition 
path sampling (TPS) to calculate r for water around Na+ 
and Fe^^; this technique is more general than the reactive 
flux method because it does not require a pre-defined 
reaction coordinate. 

By itself, the rate constant gives no information about 
the water exchange mechanism. It can, however, be used 
to obtain the activation volume (the difference in volume 
between reactant and transition state), which can be used 
to shed more light on the process. This is defined as 


( 2 ) 
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where /? = l/k^T is the inverse temperature and p the 
pressure. AV^ < 0 suggests an associative mechanism 
(with an intermediate that has a higher coordination 
number than the reactant or product), and AV^ > 0 a 
dissociative mechanis m (w hose intermediate has a lower 
coordination number) The activation volume c an be 
probed experimentally using variable-pressure NMR,ESIII1 
although this technique has not been applied to solutions 
of monovalent ions, whose water exchange rates are too 
fast for water in the hydration shell to be distinguished 
from that in the bulk. It has been used to find activation 
volumes for exchange around di- and trivale nt me tal ions, 
including transition metals and lanthanides 

Direct simulation can be used to calculate activation 
volumes for water exchange: Hermansson et al. extracted 
AV^ for Li+ and Na"*" ions from the partial mola r volu mes 
of the species involved in the exchange reactions.^^^^^Rus- 
tad and Stack computed the rate constant k at multiple 
pressures using the reactive flux method and obtained 
ARt for Li+ by differentiation.l^^ Thi s metho d has also 
been applied by several other ant hors Trajectory 

analysis of exchange events was used by Spangberg et al. 
to deduce an associative exchange mechanism around Li“'", 
which agrees with the results of other studies, without 
the need for explicit calculation of AV^^ 

Nuclear quantum effects (NQEs), such as zero-point 
energy and tunnelling, have been shown to affect the struc¬ 
tural and dynamical properties of aqueous systems 
and a natural question is to what extent this is also true 
for solvation dynamics. Videla et al. studied how quan¬ 
tum fluctuations affect the solvent relaxation following a 
sudden jump in the char ge o f a neutral solute, as in an 
electron or proton transfer. HS They found that in quantum 
mechanical simulations, the water relaxed more rapidly. 
The quantum effect on MRTs around several halide ions 
was recently investigated by Habershon, who found that 
the classical t was up to 40% larger than the quantum r, 
the ratio decreasing with increasing ion size.^^ However, 
this paper used the rigid-body SPC/E water model, and it 
has been established that neglecting flexibility and anhar- 
monicity in the 0~H bond can lead to an o verestimation 
of the significance of quantum effects.S^^^ A rigid-body 
water model allows for quantum effects in the intermolec- 
ular motion, but it does not allow for quantum effects in 
the intramolecular motion, which tend to operate in the 
opposite direction. 

A clear example of these competi^ effects is provided 
by the self-diffusion in liquid water.SS On one hand, zero 
point energy in the intermolecular modes of the liquid 
assists the breaking of hydrogen bonds and speeds up 
the diffusion. On the other, zero-point energy in the (an- 
harmonic) 0-H stretch increases both the average O-H 
bond length and the molecular dipole moment, strength¬ 
ening the hydrogen bonds and slowing down the diffusion. 
The net increase in the diffusion coefficient on including 
NQEs is therefore smaller for a flexible water model than 
a rigid-body water model.l^ 

The results of the competition between intermolecular 


and intramolecular NQEs can be subtle: Li et al. have 
found that hydrogen bonds weaker than those in pure 
water are further weakened by NQEs, while stronger 
hydrogen bonds are strengthened.^^ In terms of water 
exchange, this leads to the interesting possibility that, 
if the water-fluoride hydrogen bond is strong enough, 
quantum effects will strengthen it, making it more difficult 
for water to leave the first hydration shell of the anion. 
This would give a quantum mechanical residence time 
larger than the classical one. 

In this paper, we investigate this possibility, and also 
explore the quantum effect on the MRT in the first coordi¬ 
nation shell of a cation, by looking at the water exchange 
dynamics around Li'*' and F“ using classical molecular dy- 
namics and ring polymer molecular dynamics (RPMD).I^ 
The reason for focussing on these two ions is that they 
are the lightest in their groups and they also have the 
strongest interactions with their hydration shells: they 
are therefore the ions for which NQEs are expected to be 
largest among the alkalis and halides. Having found the 
quantum effect on t, we also look at the mechanism of the 
exchange in the two cases, and investigate the quantum 
effect on the mechanism by analyzing trajectories. 

The remainder of this paper is set out as follows: Sec- 
tion|TT]describes the models we have used for the ion-water 
interactions, as well as the various methods we have used 
for calculating the rate and mechanism of water exchange; 
Section iHll describes the results of our simulations and dis¬ 
cusses the various quantum effects observed; and Section 


IV concludes the paper. 


II. COMPUTATIONAL DETAILS 
A. Ion-Water Interaction Potentials 

We used the flexible q-TIP4P/F potential to describe 
the interactions between water molecules.^ This model is 
specifically designed for use in path integral simulations, 
in which it successfully captures several thermodynamic, 
structural, and dynamical properties of liquid water, in¬ 
cluding its density, 0-0, O-H and H-H radial distribution 
functions (RDFs), and self-diffusion coefficient The po¬ 
larizability of the water molecules is not included explicitly 
in the q-TIP4P/F model, and we chose to treat the ion- 
water interactions in the same way. As small ions, Li+ 
and F“ have small polarizabilities ,1^ so a non-polarizable 
model which takes polarization effects into account in a 
mean-field manner should be sufficient, especially since 
our goal here is simply to assess the importance of NQEs 
in the hydration of the two ions. These effects are rather 
small, so the use of a simple potential energy function that 
permits the collection of good statistics is of paramount 
importance. A polarizable water model would almost 
certainly be better for making quantitative comparisons 
with experiments,!^ but that is not our aim here. 

There are many non-polarizable ion-water interaction 
potentials available in the literature, and several different 
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TABLE I. Parameters in the present ion-water interaction 
potentials. 


Ion 

£10 [kcal/mol] 

(Tio [A] 

Potential 

Li+ 

0.136 

2.297 

Flj 


0.054 

3.791 

14-9 


sets of experimental target data to which they have been 
parameterized.^^^HMl Two popular choices of target data 
are the standard free energy of hydration, and 

the position r* of the first maximum in the ion-oxygen 
RDF.lSiHlll 

We have used these two pieces of experimental 
data to parameterize our models, and shall describe how 
they were calculated for each ion in Sec. II.B. 

The ion-water potentials consist of long-range Coulomb 
interactions between the ion and the partial charge sites 
on the water molecules, as well as short-range (van der 
Waals attraction and Pauli repulsion) interactions between 
the ion and the oxygen atoms. For Li"*", the latter were 
modelled by a Lennard-Jones potential, 


Vu{r) = 4e, 




with r the ion-oxygen distance. The parameters £10 and 
(Tio were based on those of Peng et 0 !^ with standard 
Lorentz-Berthelot combination rulea^ used to convert 
the ion-ion interaction parameters in Ref. [55] into ion- 
water parameters for the q-TIP4P/F water model. The 
lithium-oxygen potential thus obtained was found to agree 
reasonably well with both pieces of experimental data, 
without requiring any alteration. 

On the other hand, for the F“ ion, we were unable 
to parameterize the Lennard-Jones potential in such a 
way that both pieces of data were matched satisfactorily. 
A softer potential, with dependence instead of 
gave much better agreement with the target data, 


Ve,.Q{r) 



where, as in Eqn. ([^, £10 is the well depth. To parameter¬ 
ize this model, we began with the Lennard-Jones parame¬ 
ters of Jensen and Jorgensen.l^ Starting from the same 
£10 and a value of CTio that gave the potential minimum 
at the same ion-oxygen separation as in their potential, 
both values were varied to give the best agreement with 
AG^d r*. Table 1^ summarizes the resulting ion-water 
interaction parameters for both ions. 


B. Validation Calculations 

We first calculated various thermodynamic and struc¬ 
tural properties of the ionic solutions of Li+ and F“ at 


infinite dilution, and compared them with experimental 
values to validate our ion-oxygen potentials. Here, and in 
the remainder of this paper, both classical and 32-bead 
path integral simulations were carried out with 215 water 
molecules and 1 ion, at the experimental density of liquid 
water. Where temperature was held constant, it was fixed 
at 298 K using a Langevin thermostatting scheme.!^ The 
efficiency of the path integral simulations was enhanced 
by using ring-polymer contraction for the long-range part 
of the Coulomb potential.!^ 

In order to calculate r*, we ran classical molecular 
dynamics (MD) and path integral molecular dynamic^^ 
(PIMD) simulations in the canonical ensemble. For the 
Li+ ion, 25 MD simulations of length 2 ns and 25 PIMD 
simulations of length 2.6 ns were used. For both quantum 
and classical F“, 20 simulations of length 250 ps were 
carried out. From these, gio{r)^ the ion-oxygen RDF, and 
giH^r), the ion-hydrogen RDF, were extracted. The first 
minimum, of gioif) was then used to calculate the 
coordination number, 


^coord = 47rp / r‘^g,o{r)dr, (5) 

Jo 

with p the oxygen number density. While was not 

explicitly used as a target for parameterization, it did 
help us to appraise how well our models reproduced the 
experimental hydration structures of the two ions, as 
well as the extent to which quantum effects altered these 
structures. 

A raw free energy of hydration, AGhyd, was calculated 
using a thermodynamic integration scheme in which the 
hydration of the ion was broken down into two stages: 
in the first, an uncharged atom was “grown” in a box of 
215 water molecules, and in the second, this atom was 
charged.^ Throughout, the total ion-water interaction 
could be written a^^ 


Z=1 


1=1 


gigj 

47rerij ’ 


( 6 ) 

where the first sum is over all oxygen atoms, and the 
second over the charge sites on all water molecules. Here 
qi is the charge of the ion, qj the charge of site J, and rij 
the distance between the ion and site j. Ai determines 
the strength of the short-range ion-oxygen potential and 
A 2 the strength of the Coulomb interaction, so that Ai = 
A 2 = 0 corresponds to a simulation box without the ion, 
and Ai = A 2 = 1 to a fully hydrated ion. 

During the first stage of the integration, A 2 = 0 and 
Ai increases from 0 to 1, and during the second stage, 
Ai = 1 and A 2 increases from 0 to 1. A 12-point Gaussian 
quadrature was used to calculate the free energy change 
for each stage, ^ giving 24 pairs of values (Ai, A 2 ); for 
each of these, the system was evolved for 200 ps in the 
NPT ensemble. The pressure was held constant using the 
barostatting scheme of Ceriotti et a/. ,1^ which extends to 
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path integral simulations the scheme proposed earlier by 
Bussi et The thermodynamic integration for each 
ion was repeated 5 times to obtain an average value for 
and standard error in AGhyd- 

To obtain AG^^ from this raw AGhyd, we added two 
correction terms. The first accounts for the fact that 
AG^jj is the free energy difference between ions in the 

gas phase at standard pressure p®" = 1 bar and in the 
solution phase at standard molality 6”®"= 1 mol kg“^. 
This correction, equal to .RTln = 1.90 

kcal/mol, is the free energy change when an ideal gas 
at standard pressure is compressed to the molar volume 
of the ionic solution.!^ The second correction, equal to 
Nj^qicj), accounts for the p otenti al difference (j) on crossing 
the liquid-vapour interface.^^^^^ Repeating the analysis in 
Ref. I^with our water model, we found a surface potential 
of —437 mV at 298 K. This gives a correction of ±10.1 
kcal/mol to AGhyd, with the positive sign corresponding 
to F“ and the negative sign to Li+. 


sampling of W{r) around the transition state. A large 
increase in simulation time would be required to overcome 
this poor sampling, so we opted instead to use thermo¬ 
dynamic i ntegra tion to calculate W (r) in the region of 
the barrier .1221691 classical case, we ran simulations 

with the distance between the Li+ ion and the centre of 
mass of a water molecule constrained between 2.2 A and 
3.8 A, in increments of 0.025 A, and in the path integral 
case, the same constraints were applied to the centroid 
distance r. Constraints were enforced using a method 
based on RATTLE, but with the single Lagrange mul¬ 
tiplier calculated analytically.!^ The mean force on the 
constraint was recorded for each distance, and integrated 
to give a smooth curve for the PMF in the region of rC 
Transition state theory neglects dynamical recrossing 
events, which are generally quite significant in water ex¬ 
change reactions when t he ion-wa ter distance is used 
as the reaction coordinate.l2E2E3I^ These events can be 
taken into account by computing the transmission coeffi¬ 
cient K,{t)\ 


C. Rate Theory Calculations 


k = lim (9) 

t—^OO 


We followed the procedure outlined by Rey and Hynes 
to calculate the rate coefficient k for the first-order process 
of water exchange.!^ Within the transition state theory 
approximation, this rate coefficient is given by 


^TST 




(7) 


where r is the reaction coordinate (the distance between 
the ion and the centre-of-mass of the exchanging water 
molecule), the transition state separation between the 
ion-water pair, p the reduced mass of this pair, and W(r) 
the potential of mean force (PMF) along this coordinate. 
The PMF can be calculated from (?i_h 2 o(’'), the ion-water 
centre-of-mass radial distribution function, as 


W{r) = -f3 (8) 

The static simulations we used to find gio{r) and gini^) 
were also used to find 5 i_h 2 o(?'), but with one important 
difference. The correct PIMD expression for a radial dis¬ 
tribution function such as gmif) involves an average over 
the ring polymer beads, whereas the reaction coordinate 
in RPMD rate theorjEIl is more convenientljl^ chosen to 
be a function of the centroid of the ring polymer .1^21111 For 
the quantum-mechanical PMF we therefore calculated 
5 i-h2o( 7), where r is the distance between the centroids 
of the ion and the water molecule. 

While there were no problems with these calculations 
for the F“ anion (either classically or quantum mechani¬ 
cally), ions such as Li+ for which gi_H 2 o(?') is very small 
in the region of are well-known to suffer from poor 


Classically, we used the method of react ive flux to cal¬ 
culate the transmission coefficientand quantum 
mecha nical ly, ring polymer molecular dynamics rate 
theory.l^SlSfilThe RPMD prescription for /t(t) for an n-bead 
ring polymer, which reduc es to the classical prescription 

when n = 1, is simpljlSZlSS) 


n{t) = 


(r(0)5 (r(0) — r^) 9 [r{t) — r^-] 
(r(6)(5(r(0)-rt)6»[f(0)])„ 


( 10 ) 


where 6 [a;] is the Heaviside step function and (• • • )ra de¬ 
notes a Boltzmann average over the ring polymer phase 
space at reciprocal temperature Pn = P/n. 

Starting configurations for our calculations of K{t) were 
generated by running constrained simulations in the NVT 
ensemble. In the classical case, the reaction coordinate r 
was fixed at , and in the path integral case f was fixed at 
rb 15 constrained runs of 8 ns each were carried out, with 
configurations stored every 4 ps. For each run, the stored 
configurations were used to start 2000 microcanonical 
trajectories, along which the transmission coefficient was 
computed using Eqn. (10). We averaged over all 15 runs 


to find the standard error in the computed rate constant 
k, and hence in the mean residence time r = k~^. 


D. Direct Simulations 

Finally, in order to understand the mechanism of water 
exchange around Li“'" and F“ more fully, we applied the 
direct simulation method to both systems: classical calcu¬ 
lations were run in the NVE ensemble, and path integral 
simulations with the recently-developed thermostatted 
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RPMD (TRPMD) methodP In each case, we ran 20 sim¬ 
ulations of 4 ns each, and recorded every time an exchange 
event occurred. We defined an exchange to happen if a 
water molecule was in the first hydration shell of the ion 
for 5 ps in the past (except for transient escapes lasting 
no more than 2 ps), and then spent 95% of the following 
5 ps out of the shell. This meant that only very brief 
excursions back into the first shell were allowed. For each 
water molecule that was identified as having left the first 
shell, we also found the molecule that had entered to 
replace it, by looking backwards in time and using the 
same criteria to identify the molecule that was not in the 
shell before the exchange event. 

For every exchange event, the distances of the departing 
and arriving water molecules from the ion were recorded 
as a function of time, and these distances were averaged 
over events to give the (microcanonical) ensemble average 
behaviour of the exchange. From this we could infer the 
associative or dissociative nature of the reaction. 

These simulations were also used to calculate the mean 
residence time in the first coordination shell, from a resi¬ 
dence correlation function.l^ For N water molecules, this 
correlation function was defined in the quantum mechani¬ 
cal calculation as 




/ n 


( 11 ) 


Here Si{t; t*) = 1 if the centroid of the i**' water molecule 
remains in the first coordination shell of the ion between 
time 0 and time t, except for excursions out of the shell 
of duration no greater than t*. This is chosen to be 
2 ps, roughly the time that a water molecule takes to 
diffuse over one molecular radius .1^ The classical residence 
correlation function was defined in exactly the same way, 
but with n = 1 (i.e., with the centroid of the ring polymer 
of each atom in the path integral simulation being replaced 
by the position of the atom in the classical simulation). 
In both cases (quantum and classical), the MRT r was 
calculated by fitting the long-time decay of C„^(t) to 
g-t/r^ The MRT thus computed was then compared with 
the results of our reactive flux calculations. 


III. RESULTS AND DISCUSSION 
A. Validation of Potential Models 


Table [IT] compares various properties of our ion-water 
potentials, obtained from classical and path integral cal¬ 
culations, with experimental values. For the Li+ ion, the 
agreement between the results of simulations and those 
of experiment is reasonably good. Any better agreement 
is hampered by two obstacles: firstly, as pointed out by 
Jensen and Jorgensen, r* is positively correlated with both 


tJio and £ 10 , and 




is negatively correlated with 


TABLE II. Comparison of thermodynamic and structural 
properties of our ion-water models, from classical and path- 
integral simulations, with the results of experiment; the stan¬ 
dard errors in the final digits are shown in parentheses. 
is the standard free energy of hydration, r* is the first maxi¬ 
mum in the ion-oxygen radial distribution function, and ricoord 
is the coordination number of the ion. 


Property 

Classical PIMD Experiment 

Li'*' Ion 

-AG^d [kcal/mol] 
r* [A] 

^coord 

121.2(1) 120.0(1) 126.5 

1.92 1.93 2. Of 

4.03 4.05 ^ 

(1 

3 

tl 

F Ion 

-AG^d [kcal/mol] 
r* [A] 

^coord 

108.4(1) 109.2(1) 102.5 

2.63 2.63 2.61 

5.67 5.68 4-f 

(1 

B 

d] 



Ref. El 
Ref. EH 
Ref. Ei 
Ref. El 


both, precluding the simultaneous fitting of both pieces 
of data.l^ Secondly, we found that altering CTiq and Eio to 
favour one piece of data caused dramatic changes in the co¬ 
ordination number. We therefore used the Lennard-Jones 
parameters obtained from Ref. 1551 without adjustment to 
describe the solvation of lithium ions. 

For the F“ ion, we were unable to obtain a satisfactory 
fit to and r* for the q-TIP4P/F water model using 

a standard Lennard-Jones potential for the ion-oxygen 
interaction. But by varying the hard-wall part of the 
potential, and arriving finally at the dependence 
in Eqn. 0. we were able to match the experimental 
value of r* perfectly while obtaining reasonable agreement 
(as good as we were able to obtain for Li+) with the 
experimental AG^^. This is our justification for using 
the somewhat unusual 6-9 potential for the interaction 
between the fluoride ion and the water oxygens.^ 

Although the differences between the classical and quan¬ 
tum results for r* and in Table [T^ are negligible, 

nuclear quantum fluctuations do in fact have a slight ef¬ 
fect on the structure of the first hydration shell. Fig. [J 
shows the ion-oxygen and ion-hydrogen RDFs for Lu 
obtained from classical MD and PIMD simulations, and 
Fig-i shows those for F . The structures of the hy¬ 
dration shells are further illustrated by the snapshots in 
Fig.i taken from our classical MD simulations. The first 
solvation shell is arranged tetrahedrally around Li'*' and 
octahedrally around F“, with four water oxygens pointing 
towards Li+ and six water molecules donating hydrogen 
bonds to F“. For both ions, the quantum gio{r) is very 
similar to the classical one, but with a slightly smaller 
and broader first peak. The same is true for the peaks in 
ginif) corresponding to the H atoms of first-shell water 
molecules: for F“ this is both peaks, and for Li'*', the first 
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FIG. 1. Ion-oxygen (top panel) and ion-hydrogen (bottom 
panel) radial distribution functions for Li"'"(aq) from PIMD 
(solid black lines) and classical MD (dashed red lines) simula¬ 
tions at 298 K. 


FIG. 2. Ion-oxygen (top panel) and ion-hydrogen (bottom 
panel) radial distribution functions for F“(aq) from PIMD 
(solid black lines) and classical MD (dashed red lines) simula¬ 
tions at 298 K. 


one. While the quantum effect on 5 ih(^) is significantly 
larger than that on gio (r) for F“, the effects on the two 
RDFs are comparable for Li’*'. 

The decrease in intermolecular structure in quantum 
simulations of aqueous systems is well-known, and re¬ 
flects a weakening of the hydrogen- bond ing network due 
to intermolecular zero-point energy.^^^lZSl xhe ion-oxygen 
RDFs are most relevant to water exchange, and one can 
see from these that NQEs lead to a slightly less tightly 
bound first coordination shell. The fact that this effect 
is more pronounced for Li“'" than for F“ is likely due 
to a combination of two other factors. The first is the 
mitigation of water-fluoride hydrogen bond weakening by 
zero-point energy in the 0-H bond stretch, which tends 
to strengthen the hydrogen bond as discussed in Sec. I. 
Secondly, the first peak in (?Li+o(^) is sharper than that 
in gF-o(^)- This means that the water molecules are 
more tightly bound in the first shell of lithium, with more 
zero-point energy in the coordinates involved in water 
exchange.l^ We will investigate these two factors further 
in Sec. IIIIDI 

Turning now to the free energy of hydration, we note 
that while the difference between our computed quantum 
and classical values is on the order of 1-2 U-qT (which is 
larger than the uncertainty in these values), it is around 


1% of 




in each case, representing a very small 


quantum effect. Overall, quantizing these systems has 
quite a small effect on their static properties. 


B. Rate Theory Results 

In Fig. 1^ we compare the potentials of mean force for 
escape of water from the first solvation shell, obtained 
from classical and quantum simulations of both ions. The 
barrier to exchange of water around Li’*' in our classical 
simulation is 6.54 feeT, w hich is w ithin the range of 5-7 
k^T given in the literature. h^F^ l ^^ lFor F“, the correspond¬ 
ing barrier is 3.43 k^T. The heights of both barriers are 
very slightly smaller in the quantum PMFs, by 0.26 A:bT 
for Li+ and by 0.08 /cbT for F“. This decrease in the 
barrier height for F“ is less than half the decrease due to 
NQEs that HabershoiP^ found in his calculations with a 
rigid-body water model (0.2 /cbT), and we shall attribute 
this below to the absence of competing quantum effects 
in these rigid-body simulations 

The transmission coefficients K(t) obtained from our 
classical MD and RPMD simulations are shown in Fig. 

In all cases, the plateau value k is quite small. This im¬ 
plies significant recrossing of the transition state dividing 
surface, indicating that the ion-water distance r is not the 
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FIG. 3. Snapshots of the first hydration shells of Li'*' (top) 
and F“ (bottom), from classical MD simulations at 298 K. 


FIG. 4. Potentials of mean force for escape of water from the 
first solvation shells of Li^ and F“ from PIMD (solid black 
lines) and classical MD (dashed red lines) simulations at 298 

K. 


optimum reaction coordinate for these water exchange 
reactions. It is likely that a better coordinate would 
be a collective function of the positions of all the water 
molecules involved in hydration, whose rearrangement is 
necessary for a molecule to leave the first shell. 

Within the computed error bars (which are on the order 
of the widths of the lines in Fig. 5), the degree of recrossing 
in the quantum mechanical and classical simulations is 
the same. Almost all of the nuclear quantum effects on 
the rate constant are thus captured in the transition-state 
theory treatment of the reaction, and can be attributed to 


the weakened hydrogen bonding described in Sec. Ill A 


The transition state theory rate constants for 

wat er exchange are given in the first column of Table 
The quantum rate constant k^^ for each ion is 


III 


larger than the classical rate constant For Li+, 

= 1.21, and for F“ this ratio is 1.09. These 
results are entirely as predicted in Sec. Ill A the water is 
less strongly bound to the ions in the quantum simulations, 
the quantum effect being greater for Li+ than for F“. 


A similar observation was made by Habershon, who 
calculated potentials of mean force for exchange around 
halide ions, and found that a smaller ion-water binding 
strength was correlated with a smaller decrease in barrier 
height due to zero-point energy.!^ Although the Li'*' and 
F“ ions interact differently with the molecules in their 
first hydration shells, they appear to fit into this pattern. 


Since water will be more tightly bound to Li"^ than to the 
other alkali metals, and the other alkalis are also heavier 
than lithium, we would predict that the quantum effects 
on seen here are the largest that will be seen in this 
group. 

The quantum and classical transmission coefficients 
K, recrossing-corrected rate coefficients k = and 

mean residence times r = k~^, are also given in Table 
m According to our RPMD calculations, the MKT of 
water in the first hydration shell of Li+ is 118 ps, and in 
the first shell of F“, 31 ps. There are no experimental 
results with which we can directly compare these values, 
but neutron scattering experiments suggest that t < 100 
ps for both ionsFor F“, our results are comfortably 
within this bound, while for Li+ the imprecision in the 
experimental result makes our value reasonable. The 
MKTs of Li'*' from earlier computational studies span 
two ord ers of magnit ude, with values from 25-400 ps 
reported i ^b2 | 23 | 24 | 29 | 79 | those of F are in the 17-25 

ps range P^SHHini 


C. Direct Simulation Results 

The mean residence times calculated by fitting the 
residence correlation function to an exponential 














TABLE III. Dynamical properties of water exchange around Li"*" and F“ from classical and quantum simulations, is the 

transition state theory rate constant, k is the transmission coefficient and k = is the full rate constant. Trf = 1/fc is the 

MRT calculated using the reactive flux method and Tds is the decay time of the residence correlation function calculated from 
direct simulations, rd/rqm is the ratio of classical to quantum mean residence times. The standard errors in the final digits are 
shown in parentheses. 




K 

k [ns 

Trf [ps] 

Tds [ps] 

Li+ Ion 

Classical 

42 

0.15(1) 

6.2(4) 

160(11) 

160(9) 

Quantum 

51 

0.16(1) 

8.2(5) 

121(8) 

116(4) 




'Tcl/ '7”qm 

1.34 

1.38 

F Ion 

Classical 

292 

0.10(1) 

29.2(2.9) 

34.2(3.4) 

29(1) 

Quantum 

319 

0.10(1) 

31.9(3.2) 

31.3(3.1) 

26(1) 




"^cl/ "^qm 

1.09 

1.11 




FIG. 5. Transmission coefficients for the water exchange 
reactions around Li'*' and F“ from RPMD (solid lines) and 
classical MD (dashed red lines) simulations. 


decay are given in the final column of Table |III[ They 
compare reasonably well with the results of our rate theory 
calculations, which gives us confidence in our calculated 
values. Similar agreement between the two methods was 
also found by Spangberg et al^ and by Kerisit and 
Rosso.ti^ 

To complement our calculations of r, Fig. shows the 
average ion-water distance for the leaving water molecule 


and for its replacement during an exchange event. From 
this, we can classify the mechanism of exchange around 
the ion. According to Langford and Gray,l^ ligand ex¬ 
changes may be associative, A (in which the initial step 
is the attachment of the replacing ligand), dissociative, 
D (in which the initial step is the detachment of a lig¬ 
and), or interchange, I (with the ligands exchanging in a 
concerted fashion). This latter mechanism can be further 
split up depending on whether it has more associative or 
dissociative character (R and Id respectively). 

Using these classifications, we look first at exchange 
around Li'*': the leaving and replacing water molecules 
cross within the first solvation shell, but the leaving 
molecule has begun its exit from the shell by the time they 
cross. This corresponds to an associative interchange (U) 
mechanism, which agrees with the results of Spangberg et 
al.,l^who found that the majority of exchanges around 
Li'*' had an associative nature (either A or U). The acti¬ 
vation volumes calculated in the literature are n egative, 
which also points to an A or U mechanism.E^Ell 

For F“, averaging over all exchange events suggests that 
the two molecules cross outside the shell, but very close to 
its boundary, in a dissociative interchange (Id) mechanism. 
However, Heuft and Meijer found that exchange around 
F“ proceeds via two pathways, one of which is associative 
and the other dissociative.^^ To investigate whether the 
same was true for our model, and the behaviour in Fig. 
was thus an average of two mechanisms, we took each 
individual event and classified it as associative or disso¬ 
ciative, based on whether the molecules crossed inside or 
outside of the first shell. 

We found that about half of the events were associative, 
and half dissociative, suggesting that there are indeed 
two pathways to exchange. In Fig. we show the dis¬ 
tances of the entering and leaving molecules from F“, 
separately averaged over the associative and the disso¬ 
ciative exchanges. The distinction between the two is 
clear-cut: in both mechanisms, which are classified as U 
and Id respectively, the molecules cross further from the 





















9 



FIG. 6. Average trajectories of entering and leaving water 
molecules during exchange events around Li"*" and F“ from 
TRPMD (solid black lines) and classical MD (dashed red 
lines) simulations. In each case, the dotted black line shows 
the position of the first shell boundary. 


first shell boundary than in Fig.|^ Repeating the analysis 
for Li+, we found that fewer than 3% of the exchanges 
had any dissociative character, so exchange around this 
ion can safely be described as associative. 

In both Figs. and the results of our TRPMD simu¬ 
lations are compared to those of classical MD. It is evident 
that, once again, the difference between the quantum and 
classical simulations is minimal. Coupled with our results 
from previous sections, we see that the overall effect of 
quantum fluctuations on water exchange around these 
two ions is very small. Quantum fluctuations do slightly 
enhance the rate of water exchange in both cases, by 
^ 10% for F“ and ~ 35% for Li+ (see Table III), but 
this enhancement can be traced almost entirely to the 
reduced barrier heights of the quantum mechanical poten¬ 
tials of mean force in Fig. 4. Quantum effects have very 
little impact on either the time-dependent transmission 
coefficients in Fig. 5 or the average reactive trajectories 
in Figs. 6 and 7. Presumably this is because these are 
determined by the dynamics at configurations on either 
side of the exchange barrier, where there is hardly any 
difference between the quantum and classical potentials 
of mean force (see Fig. 4). 



FIG. 7. Average trajectories of entering and leaving water 
molecules during associative and dissociative exchange events 
around F~ from TRPMD (solid black lines) and classical MD 
(dashed red lines) simulations. In each case, the dotted black 
line shows the position of the first shell boundary. 


D. Competing Quantum Effects 

There is still one remaining aspect of our results that 
needs explaining. We have shown that nuclear quantum 
effects decrease the MRT of water in the first coordina¬ 
tion shell of F“ by only ^ 10%, which is considerably 
smaller than the 40% decrease found by Habershon.^^ We 
shall now establish that this is because of the competing 
quantum effects that are present in our simulations: the 
0-H- • • F“ hydrogen bond is strengthened by zero-point 
energy in the anharmonic O-H stretch, and weakened 
by zero-point energy in perpendicular directions. The 
first of these effects is missing from the rigid-body water 
simulations of Ref. [28l 

A convenient probe of competing quantum effects is 
provided by the quantum kinetic energy (T). If this is 
greater for the H atoms hydrogen-bonded to F“ than the 
H atoms in bulk water, the former will be more confined 
and have a higher zero-point energy than those in the 
bulk.l^ What makes this especially useful is that (T) 
can be resolved into contributions from motion in three 
orthogonal directions, which provides a very natural way 
to separate the competing quantum effects.!^ 

As in Ref. HU we calculated the centroid virial 
estimatoi!^ for the kinetic energy tensor, Tap = PaPp/^m, 
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TABLE IV. Components of the quantum kinetic energy tensor for H atoms in various environments relative to the Li'^ and 
F“ ions in PIMD simulations at 298 K. The left-hand panel illustrates the tensor, with the directions of its three eigenvectors 
shown, and (T) is the sum of the three eigenvalues. 


-H 



[meV] 

(Ti) 

(T 2 ) 

(Ta) 

(T) 

Li+ Ion 

First shell water 

98.0 

34.3 

21.8 

154.1 

Bulk water 

98.8 

33.7 

21.6 

154.1 

F Ion 

0-H- • • F“ hydrogen bonding 

91.6 

38.5 

28.2 

158.3 

Other first shell hydrogen 

99.0 

33.8 

21.6 

154.4 

Bulk water 

98.9 

33.7 

21.5 

154.1 


for each H atom in each snapshot of our PIMD simulation. 
We then found the rotation that minimized the mean- 
square difference between the water molecule containing 
the H atom and a reference molecule, and applied this 
rotation to the kinetic energy tensor, so as to obtain the 
tensor in the molecular frame. The resulting tensors were 
averaged over the simulation for three different groups of H 
atoms: those in the first hydration shell hydrogen-bonded 
to F“, those in the first hydration shell not hydrogen- 
bonded to F“, and those in the remainder of the solution 
(described here as bulk). For each of these species, the 
averaged tensor (Tc/j) was then diagonalised to obtain its 
eigenvalues (T^) and the corresponding eigenvectors. 


Table IV reports the total kinetic energies (T) and 
their components (Ti) for each of the three types of H 
atom. In all cases, the components are considerably larger 
than the classical expectation value of k-QT/2 = 12.8 
meV, because of the significant amount of zero-point 
energy in 0-H bonds. One sees that the total kinetic 
energy (T) differs appreciably from its bulk value only for 
those atoms directly hydrogen-bonded to the fluoride ion. 
(Ti), the component in the 0-H direction, has decreased 
from its value in bulk water, indicating a decrease in 
confinement in this direction due to stretching of the 
bond to facilitate hydrogen bonding. On the other hand, 
(T 2 ) and (T 3 ), the components perpendicular to the O- 
H bond, have increased, reflecting more confinement in 
these directions due to the stronger hydrogen-bonding 
interaction of H with the F“ ion than with the oxygen 
atoms of other water molecules. This is a clear indication 
of competing quantum effects. The net result is that (T) 
is larger by 4.2 meV^ U-qT/Q for the H atoms hydrogen 
bonded to F“ than the H atoms in bulk water. The first 
solvation shell of the F“ is therefore destabilised slightly 
by quantum effects, which is consistent with the slight (~ 
10%) reduction in the MRT of the water molecules in this 
shell we have found in our simulations. Without allowing 
for the anharmonicity of the 0-H bond, the decrease 
in the (Ti) component between the bulk and the first 
solvation shell would be expected to be smaller, resulting 
in a larger increase in both (T) and the exchange rate on 
including NQEs, as Habershon found in his simulations .1^ 


Table IV also compares the components of (T) for H 
atoms in the first shell of Li+ with those in the bulk. The 
total kinetic energy is the same in both environments, 
though there is a small difference in the individual com¬ 
ponents: as for F“, (Ti) is smaller in the first shell, while 
(T 2 ) and (T 3 ) are larger. However, in this case the kinetic 
energy of the H atoms does not play any part in explain¬ 
ing the quantum effect on the exchange rate, since the H 
atoms do not participate in hydrogen bonding with the 
ion. Our results for Li+ are in agreement with the results 
of Videla et a/.,S^who found that the spatial delocaliza¬ 
tion of H atoms in the first solvation shells of cations was 
very similar to that in bulk water. From the point of view 
of the H atom kinetic energy, there is certainly very little 
distinction between the two environments. 


The correlation between halide-water binding strengths 
and quantum effect on MRTs has been attributed 
to greater zero-point energy in more strongly bound 
systems.l^ Since Li+ and F“ interact differently with 
their first hydration shells, it is not immediately clear 
whether or not this factor contributes to the larger quan¬ 
tum effect for the former. For this reason, we have also 
used transition state theory to look at water exchange 
around water. This molecule interacts with its first shell 
as both a hydrogen bond donor and acceptor, so is in this 
sense intermediate between the two ions. Computing the 
PMF for this exchange using MD and PIMD simulations 
with 216 water molecules, we found a barrier height of 
1.4 k^T, even smaller than that for F“ (3.4 k^T). In a 
similar vein, is equal to 1.61 ps“^ classically and 

1.69 ps“^ quantum mechanically: thi s incre ase of 5% is 
also smaller than the 9% found in Sec. IlHBl for F“. This 
lends support to the suggestion that, regardless of the 
mode of solvent-solute binding, more strongly bound sol¬ 
vent molecules lead to a greater nuclear quantum effect 
on T. 


IV. CONCLUSIONS 

In this paper, we have studied the impact of nuclear 
quantum effects on the static and dynamical properties 
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of water in aqueous solutions of Li’*' and F“ ions using 
classical and path integral molecular dynamics techniques. 
We have found that for both ions, quantization causes a 
small destabilization of the first hydration shell, leading 
to a shorter mean residence time for the water molecules 
in this shell. This is explained in both cases by a net 
disruption of the hydrogen bonding on the inclusion of 
zero-point energy. The effect is greater for Li+ than for 
F“ for two reasons: Li+ interacts more strongly with the 
water molecules in its first hydration shell, and there is a 
competition between quantum effects in the 0-H- • • F“ 
hydrogen bond as discussed in Sec. III.D. Since an essen¬ 
tial ingredient in this competition is the anharmonicity of 
the 0-H bond,l^ it is missed in simulations of rigid-body 
water molecules, which explains the difference between 
our results for fluoride and those of Ref. HSj 

We have also found that the mechanism of water ex¬ 
change around both ions, exemplified by the trajectories 
of the leaving and replacing molecules, is practically un¬ 
changed by nuclear quantum effects. The real upshot of 
our results is thus that classical simulations provide a 
qualitative, and even semiquantitative, description of the 
water exchange dynamics. Since quantum effects have 


been seen to diminish with weaker ion-water binding, the 
Li'*" and F“ ions bind water the most strongly in their 
respective groups, and they are also the lightest elements 
in their groups, we would expect this to be true for all 
alkalis and halides. The classical treatment of the nuclear 
motion in most previous work on water exchange around 
these ions is therefore justified. 


ACKNOWLEDGEMENTS 

We would like to thank Anne Wilkins and Josh More for 
critical reading of the manuscript and helpful comments. 
D.M.W. acknowledges funding from the Oxford Univer¬ 
sity Clarendon Fund and St. Edmund Hall, and computer 
time on the University of Oxford Advanced Research Com¬ 
puting (ARC) facility and the IRIDIS High Performance 
Computing facility. L.X.D. acknowledges funding from 
the U.S. Department of Energy, Office of Science, Office 
of Basic Energy Sciences, Division of Chemical Sciences, 
Geosciences, and Biosciences. 


^ K. Modig, E. Liepinsh, G. Otting, and B. Halle, J. Am. 
Chem. Soc. 126, 102 (2004) 

^ D. Russo, R. K. Murarka, J. R. D. Copley, and T. Head- 
Gordon, J. Phys. Chem. B 109, 12966 (2005) 

® W. H. Casey, Chem. Rev. 106, 1 (2006) 

^ A. F. Panasci, C. A. Ohlin, S. J. Harley, and W. H. Casey, 
Inorg. Chem. 51, 6731 (2012) 

® T. K. Udeigwe, P. N. Eze, J. M. Teboh, and M. H. Stietiya, 
Environ. Int. 37, 258 (2011) 

® R. W. Impey, P. A. Madden, and I. R. McDonald, J. Phys, 
Chem. 87, 5071 (1983) 

' A. E. Merbach, Pure Appl. Chem. 59, 161 (1987) 

® H. Ohtaki and T. Radnai, Chem. Rev. 93, 1157 (1993) 

® R. Rey and J. T. Hynes, J. Phys. Chem. 100, 5611 (1996) 
L. Helm and A. E. Merbach, J. Chem. Soc., Dalton Trans, 
, 633 (2002) 

L. Helm and A. E. Merbach, Chem. Rev. 105, 1923 (2005) 
J. R. Rustad and A. G. Stack, J. Am. Chem. Soc. 128, 
14778 (2006) 

S. Kerisit and K. M. Rosso, J. Chem. Phys. 131, 114512 
(2009) 

L. X. Dang and H. V. R. Annapureddy, J. Chem. Phys. 
139, 084506 (2013) 

Z. Luz and R. G. Shulman, J. Chem. Phys. 43, 3750 (1965) 

L. Helm, Coord. Chem. Rev. 252, 2346 (2008) 

E. Balogh and W. H. Casey, Prog. Nucl. Magn. Reson. 
Spectrosc. 53, 193 (2008) 

“ P. S. Salmon, W. S. Howells, and R. Mills, J. Phys. C 20, 
5727 (1987) 

P. S. Salmon, G. J. Herdman, J. Lindgren, M. C. Read, and 

M. Sandstrom, J. Phys. Condens. Matter 1, 3459 (1989) 

P. S. Salmon, M. C. Bellissent-Funel, and G. J. Herdman, 
J. Phys. Condens. Matter 2, 4297 (1990) 


D. E. Moilanen, D. Wong, D. E. Rosenfeld, E. E. Fenn, 
and M. D. Fayer, Proc. Natl. Acad. Sci. USA 106, 375 
(2009). 

J. Minbiao, M. Odelius, and K. J. Gaffney, Science 328, 
1003 (2010). 

S. Koneshan, J. C. Rasaiah, R. M. Lynden-Bell, and S. H. 
Lee, J. Phys. Chem. B 102, 4193 (1998) 

S. H. Lee and J. C. Rasaiah, J. Phys. Chem. 100, 1420 
(1996) 

A. Luzar and D. Chandler, Phys. Rev. Lett. 76, 928 (1996) 
A. Luzar and D. Chandler, Nature 379, 55 (1996) 

A. Chandra, J. Phys. Chem. B 107, 3899 (2003) 

S. Habershon, Phys. Chem. Chem. Phys. 16, 9154 (2014) 
D. Spangberg, R. Rey, J. T. Hynes, and K. Hermansson, 
J. Phys. Chem. B 107, 4470 (2003) 

T. S. Hofer, H. T. Tran, C. F. Schwenk, and B. M. Rode, 

J. Comput. Chem. 25, 211 (2004) 

C. B. Messner, T. S. Hofer, B. R. Randolf, and B. M. Rode, 
Phys. Chem. Chem. Phys. 13, 224 (2011) 

D. Laage and J. T. Hynes, J. Phys. Chem. B 112, 7697 
(2008). 

C. H. Langford and H. B. Gray, Ligand Substitution Pro¬ 
cesses (W. A. Benjamin: New York, 1966). 

Y. Ducommun, K. E. Newman, and A. E. Merbach, Inorg, 
Chem. 19, 3696 (1980); A. D. Hugi, L. Helm, and A. E. 
Merbach, Helv. Chim. Acta 68, 508 (1985), Inorg. Chem. 
26, 1763 (1987) 

C. Cossy, L. Helm, and A. E. Merbach, Inorg. Chem. 28, 
2699 (1989) 

D. Spangberg, M. Wojcik, and K. Hermansson, Chem. 
Phys. Lett. 276, 114 (1997) 

K. Hermansson and M. Wojcik, J. Phys. Chem. B 102, 
6089 (1998) 





12 


H. V. R. Annapureddy and L. X. Dang, J. Phys. Chem. B 
118, 7886 (2014) 

A. Wallqvist and B. Berne, Chem. Phys. Lett. 117, 214 
(1985) 

T. F. Miller and D. E. Manolopoulos, J. Chem. Phys. 122, 
184503 (2005). 

F. Paesani and G. A. Voth, J. Phys. Chem. B 113, 5702 
(2009) 

J. Liu, R. S. Andino, C. M. Miller, X. Chen, D. M. Wilkins, 
M. Ceriotti, and D. E. Manolopoulos, J. Phys. Chem. C , 
2944 (2013) 

P. E. Videla, P. J. Rossky, and D. Laria, J. Chem. Phys, 
139, 164506 (2013) 

S. Habershon, T. E. Markland, and D. E. Manolopoulos, 

J. Chem. Phys. 131, 024501 (2009) 

T. E. Markland and B. J. Berne, Proc. Natl. Acad. Sci. 
USA 109, 7988 (2012) 

“ X.-Z. Li, B. Walker, and A. Michaelides, Proc. Natl. Acad, 
Sci. USA 108, 6369 (2011) 

S. Habershon, D. E. Manolopoulos, T. E. Markland, and 

T. F. Miller, Annu. Rev. Phys. Chem. 64, 387 (2013). 

G. D. Mahan, Phys. Rev. A 22, 1780 (1980); C. Hattig 
and B. A. Hefi, J. Ghem. Phys. 108, 3863 (1998) 

X. Huafeng, H. A. Stern, and B. J. Berne, J. Phys. Chem. 
B 106, 2054 (2002). 

D. E. Smith and L. X. Dang, J. Chem. Phys. 100, 3757 
(1994) 

K. P. Jensen and W. L. Jorgensen, J. Chem. Theory Corn- 
put. 2, 1499 (2006). 

I. S. Joung and T. E. Cheatham, J. Phys. Chem. B 112, 
9020 (2008) 

H. Yu, T. W. Whitfield, E. Harder, G. Lamoureux, 

I. Vorobyov, V. M. Anisimov, A. D. Mackerell, and B. Roux, 

J. Ghem. Theory Gomput. 6, 774 (2010) 

M. M. Reif and P. H. Hunenberger, J. Ghem. Phys. 134, 
144103 (2011) 

T. Peng, T.-M. Ghang, X. Sun, A. V. Nguyen, and L. X. 
Dang, J. Mol. Liq. 173, 47 (2012), note that the values of 
Vm given in Table 2 of this reference are actually rm/2 for 
the ion-ion Lennard-Jones potential. In particular, an = 
1.438 A for Li’^. 

M. P. Allen and D. J. Tildesley, Computer Simulations of 
Liquids (Clarendon Press, Oxford, 1987). 

M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. 
Manolopoulos, J. Chem. Phys. 133, 124104 (2010). 

T. E. Markland and D. E. Manolopoulos, Chem. Phys. Lett, 
464, 256 (2008) 

M. Parrinello and A. Rahman, J. Chem. Phys. 80, 860 
(1984) 


D. Horinek, S. 1. Mamatkulov, and R. R. Netz, J. Chem. 
Phys. 130, 124507 (2009) 

M. Ceriotti, J. More, and D. E. Manolopoulos, Comput. 
Phys. Commun. 185, 1019 (2013) 

G. Bussi, T. Zykova-Timan, and M. Parrinello, J. Ghem. 
Phys. 130, 074101 (2009). 

G. Lamoureux and B. Roux, J. Phys. Chem. B 110, 3308 
(2006) 

M. A. Wilson, A. Pohorille, and L. R. Pratt, J. Chem. 
Phys. 88, 3281 (1988) 

““ 1. R. Craig and D. E. Manolopoulos, J. Chem. Phys. 122, 
84106 (2005). 

I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. 123, 
34102 (2005). 

®^ R. Collepardo-Guevara, 1. R. Craig, and D. E. Manolopou¬ 
los, J. Chem. Phys. 128, 144502 (2008). 

®® T. E. Markland, S. Habershon, and D. E. Manolopoulos, 

J. Chem. Phys. 128, 194506 (2008). 

®® E. Guardia, R. Rey, and J. Padro, Ghem. Phys. 155, 187 
(1991) 

H. G. Andersen, J. Comp. Phys. 52, 24 (1983) 

Charles H. Bennett, in Algorithms for Chemical Compu¬ 
tations, ACS Symposium Series, Vol. 46, edited by R. E. 
Christoffersen (American Chemical Society, Washington, 
D. C., 1977) pp. 63-97. 

D. Chandler, J. Chem. Phys. 68, 2959 (1978) 

M. Rossi, M. Ceriotti, and D. E. Manolopoulos, J. Chem. 
Phys. 140, 234116 (2014) 

M. D. Tissandier, K. A. Cowen, W. Y. Feng, E. Gundlach, 
M. H. Cohen, A. D. Earhart, J. V. Coe, and T. R. Tuttle, 
J. Phys. Chem. A 102, 7787 (1998) 

Y. Marcus, Chem. Rev. 88, 1475 (1988) 

^® S. Varma and S. B. Rempe, Biophys. Chem. 124, 192 
(2006) 

" A few other publications have used non-polarizable ion- 
water potentials other than the Lennard-Jones. See, for 
example. Refs. l6land[29l 

R. A. Kuharski and P. J. Rossky, Chem. Phys. Lett. 103, 
357 (1984) 

G. 1. Szasz and K. Heinzinger, J. Chem. Phys. 79, 3467 
(1983) 

J. M. Heuft and E. J. Meijer, J. Chem. Phys. 122, 094501 
(2005) 

M. Ceriotti and D. E. Manolopoulos, Phys. Rev. Lett. 109, 
100604 (2012) 









